Optimal modularity for nucleation in network-organized Ising model 
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We study nucleation dynamics of Ising model in a topology that consists of two coupled random 
networks, thereby mimicking the modular structure observed in real-world networks. By introducing 
a variant of a recently developed forward flux sampling method, we efficiently calculate the rate 
and elucidate the pathway for nucleation process. It is found that as the network modularity 
becomes worse the nucleation undergoes a transition from two-step to one-step process. Interestingly, 
the nucleation rate shows a nonmonotonic dependency on the modularity, in which a maximal 
nucleation rate occurs at a moderate level of modularity. A simple mean field analysis is proposed 
to qualitatively illustrate the simulation results. 
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I. INTRODUCTION 

In the last decade, critical phenomena in complex net- 
works have received an enormous amount of attention 
in the field of statistical physics and many other disci- 
plines (see, for example, a recent review |l|). Extensive 
research interests have focused on the onset of critical 
behaviors in diverse network topology, which included a 
wide range of issues: percolation phenomenon [2j-[5j, epi- 
demic thresholds [a, Mh order-disorder transitions \t 
synchronization |13lll4j |. self-organized criticality [If 
nonequilibrium pattern formation [l7j . etc. However, 
there is much less attention paid to the dynamics/kinetics 
of phase transition itself in complex network, such as nu- 
cleation in a first-order phase transition. 

Nucleation is a fluctuation-driven process that initi- 
ates the decay of a metastable state into a more stable 
one [l8j]. A first-order phase transition usually involves 
the nucleation and growth of a new phase. Many impor- 
tant phenomena in nature, including crystallization [l9j, 
glass formation [20|, and protein folding [2l|, etc., are 
associated with nucleation. Despite its apparent impor- 
tance, many aspects of nucleation process are still un- 
clear and deserve more investigations. The Ising model, 
which is a paradigm for many phenomena in statistical 
physics, has been widely used to study the nucleation 
process. Despite its simplicity, Ising model has made im- 
portant contributions to the understanding of nucleation 
phenomena in equilibrium systems and is likely to yield 
important insights also for nonequilibrium systems. In 
two-dimensional lattices, for instance, shear can enhance 
the nucleation rate and the rate peaks at an intermedi- 
ate shear rate f22j, a single impurity may considerably 
enhance the nucleation rate (23[, and the existence of 
a pore may lead to two-stage nucleation and the overall 
nucleation rate can reach a maximum level at an interme- 
diate pore size [24| ■ Nucleation pathway of Ising model 
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in three-dimensional lattice has also been studied using 
transition path sampling approach [25]. In addition, Ising 
model has been frequently used to test the validity of 
classical nucleation theory (CNT) 2(i .'>(J. However, all 
these studies are limited to regular lattices in Euclidean 
space. Since many real systems can be properly mod- 
eled by network-organized structure, it is thus natural to 
ask how the topology of a networked system affects the 
nucleation process of Ising model. 

In a recent work [31(, we have studied nucleation dy- 
namics on scale- free networks, in which we found that nu- 
cleation starts from, on average, nodes with more lower 
degrees, and the rate for nucleation decreases exponen- 
tially with network size and the size of critical nucleus 
increase linearly with network size, implying that nucle- 
ation is relevant only for a finite-size network. Herein, 
we want to study nucleation dynamics of Ising model 
in a modular network composed of two coupled random 
networks. It is known that many real- world networks, 
as diverse as from social networks to biological networks, 
have found to exhibit modularity structures [321 . l33j , that 
is, links within modules are much more denser than those 
between modules. Many previous studies have revealed 
that such modular structures have a significant impact 
on the dynamics ta king place on the networks, such as 
synchronizatio n [33.13511 neural excitability f3(| , sp read- 
ing dynamics [37l. l3cf, opinion formation [39l l40j. and 
Ising phase transition [4lM43| . In particular, for major- 
ity model [HI and Ising model [41H43T ] in modular net- 
works, it has been shown that there exists a region in a 
discontinuous transition where modular order phase and 
global order phase coexist. However, these studies mainly 
focused on phase diagrams in parameters space, and did 
not make detailed investigation for the transition process 
from a phase to another that may undergo a nucleation 
process. 

Since nucleation is an activated process that occurs ex- 
tremely slow, brute-force simulation is prohibitively ex- 
pensive. To overcome this difficulty, we will use a variant 
of a recently developed simulation method, forward flux 
sampling (FFS) This method allows us to calculate 



2 



nucleation rate and determine the properties of ensemble 
toward nucleation pathways. We found that as the degree 
of network modularity decreases nucleation goes through 
a transition from two-step to one-step process, and the 
rate exhibits a maximum at an intermediate degree of 
modularity. Free energy profiles for different modularity 
obtained by umbrella sampling (US) [45| and a simple 
mean-held (MF) analysis help us understand the FFS 
results. 

This paper is organized as follows. In SecJIll we de- 
scribe the details of our model and the simulation method 
applied to this system. In Sec JIIIL we present the results 
for the nucleation rate and pathway in modular networks. 
In Sec lIVl a simple mean held analysis is used to qualita- 
tively illustrate the simulation results. At last, discussion 
and main conclusions are addressed in SecfVl 



II. MODEL AND SIMULATION DETAILS 

Consider a network consisting of N nodes arranged 
into two modules with N\ and N2 nodes. For simplicity, 
we only consider the case of N\ — N2 = N/ 2 throughout 
this paper. The connection probability between a pair 
of nodes belonging to the same module is pi , while that 
for nodes belonging to different modules is p a . The pa- 
rameter a = — 6 [0, 1], defined as the ratio of inter- to 
intra-modular connectivity, measures the degree of mod- 
ularity. The higher degree of modularity of a network is, 
the smaller value of a is. As a — > 0, the network be- 
comes two isolated clusters, while as a — >• 1, the network 
approaches a Erdos-Renyi (ER) random network. When 
a is varied the total number of links of the network is 
kept unchanged, N %°' , where (k) is the average degree. 

This restriction leads to pi = y7^^ and p a = j^pfe ■ 
Each node is endowed with an Ising spin variable Sj that 
can be either +1 (up) or —1 (down). The Hamiltonian 
of the system is given by 

H = -.7 f ll}" (1) 

% 

where J(> 0) is the coupling constant and h is the exter- 
nal magnetic field. The elements of the adjacency matrix 
of the network take ay = 1 if nodes i and j are connected 
and dij = otherwise. 

Our simulation is performed by Metropolis spin-flip 
dynamics [4(| . in which we attempt to flip each spin once, 
on average, during each Monte Carlo (MC) cycle. In 
each attempt, a randomly chosen spin is flipped with the 
probability min(l, e~ l3AE ), where (3 = l/(fcsT) and ks 
is the Boltzmann constant and T is the temperature, and 
AE is the energy difference due to the flipping process. 
Here, we set J = 1, h > 0, T < T c (T c is the critical 
temperature), and start with a metastable state in which 
Si = — 1 for most of the spins. The system will stay in 
that state for a significantly long time before undergoing 
a nucleation transition to a more stable state with most 



spins pointing up. We are interested in the pathway and 
rate for this nucleation process. 

FFS method has been used to calculate rate constants 
and transition paths for rare events in equilibrium and 
nonequilibrium systems [22T - I24] . HJ, H3, EH . This method 
uses a series of interfaces in phase space between the ini- 
tial and final states to force the system from the initial 
state A to the final state B in a ratchet-like manner. Be- 
fore the simulation begin, an order parameter A is first 
defined, such that the system is in state A if A < Ao and 
it is in state B if A > Am- A series of nonintersecting 
interfaces Xi (0 < i < M) lie between states A and B, 
such that any path from A to B must cross each inter- 
face without reaching Aj+i before A^. The algorithm first 
runs a long-time simulation which gives an estimate of 
the flux <&a,o escaping from the basin of A and generates 
a collection of configurations corresponding to crossings 
of interface Ao- The next step is to choose a configura- 
tion from this collection at random and use it to initiate 
a trial run which is continued until it either reaches Ai 
or returns to Ao- If Ai is reached, store the configuration 
of the end point of the trial run. Repeat this step, each 
time choosing a random starting configuration from the 
collection at Ao- The fraction of successful trial runs gives 
an estimate of of the probability of reaching Ai without 
going back into A, P(Ai|Ao). This process is repeated, 
step by step, until Am is reached, giving the probabilities 
P (Aj+i I Aj) (i = 1, • • • , M— 1). Finally, we get the transi- 
tion rate R from A to B, which is the product of the flux 
§A,o and the probability P (A A /|A ) = Hfio 1 P (<Vhi|Aj) 
of reaching Am from Ao without going into A. The de- 
tailed descriptions of FFS method see Ref . [49( . 

However, conventional FFS method will become very 
inefficient if one intermediate metastable state exists be- 
tween initial state and final state, as a two-step nucle- 
ation process demonstrated in FigjTJ This is because 
that sampling paths will be trapped in these long-lived 
metastable states so that they hardly return to initial 
state. To solve this problem, we will perform instead 
two-step samplings from initial down-spin state to inter- 
mediate metastable state, and then to final up-spin state, 
giving the two-step rates, R\ and R2, respectively. Since 
the total mean time for nucleation is simply the sum of 
the mean time of the two-step process, the total rate can 
be expressed as R = (itf 1 + -R^ 1 ) . To determine the 
location of the intermediate state, during FFS we moni- 
tor the sampling time for the probability P (Aj+i |Aj) be- 
tween two neighboring interfaces. If the sampling time 
spent on between interfaces i and i + 1 is much more 
than its previous step and the probability P(Ai+i|Aj) is 
nearly one, we consider the ith interface as the location 
of the intermediate state. If such conditions do not meet 
during the whole sampling, the intermediate metastable 
state does not exist, meaning that nucleation is a one- 
step process. Note that the method is straightforward to 
generalize to a multi-step nucleation process. 

Here, we define the order parameter A as the total num- 
ber of up spins in the networks. The spacing between in- 
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FIG. 1: (Color online) Typical time evolutions of the num- 
ber of up spins A. It is shown that the system undergoes a 
two-step nucleation process for a — 0.011 and a one-step nu- 
cleation process for a = 0.051. The representative network 
configurations at different moments indicated by arrows are 
shown in Fig[2] Other parameters are N = 400, (k) = 6, 
T = 2.0, and h = 1.2. 
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FIG. 2: (Color online) Five representative network configu- 
rations at different moments indicated in Fig[T] where down- 
spin nodes and up-spin nodes are denoted by blue circles and 
red triangles, respectively, (a)-(c) correspond to the case of 
a = 0.011 and (d)-(e) correspond to the case of a = 0.051. 



terfaces is fixed at 3 up spins, but the computed results 
do not depend on this spacing. We perform 1000 trials 
per interface for each FFS sampling, from which at least 
100 configurations are saved at each interface in order to 
investigate the statistical properties of an ensemble of re- 
active pathways to nucleation. The results are obtained 
by averaging over 10 independent FFS samplings and 50 
different network realizations. 



III. RESULTS 

To begin with, in FigUwe exhibit typical time evolu- 
tions of the number of up spins A corresponding to two 
different values of network modularity a = 0.011 and 
a = 0.051 via brute-force simulations, with relevant pa- 
rameters N = 400, (k) = 6, T = 2.0, and h = 1.2. It is 
clearly observed that the system undergoes a two-step nu- 
cleation process for a = 0.011 and a one-step nucleation 
process for a = 0.051. We also plot several represen- 
tative configurations in FigJS] corresponding to different 
phases of the system, respectively. Before the nucleation 
happening, the system lies in a metastable state, where 
most of the nodes are in down-spin state (indicated by 
blue circles), as shown in Figj2ja) and Figj2{d). When 
the network modularity is very good, the system enters 
into an intermediate metastable state via the first-step 
nucleation, where nodes in one of modules are in up- 
spin state (indicated by red triangles), while nodes in 
the other module are still in down-spin state, as shown in 
Figj2jb). When the network modularity becomes worse, 
such an intermediate metastable state disappears so that 
the nucleation becomes a one-step process. Finally, the 
system will enter into the most stable state, where al- 
most all spins are in up-spin state, as shown in Figj^Jc) 
and Figj2]Je). Moreover, we note that that the nucleation 
process typically takes the order of 10 6 or more MC steps 
that is computationally costly. Therefore, in what follows 
we will give the results obtained by FFS method. 

The nucleation rate J? as a function of a is plotted 
in Figf3fa), with relevant parameters being the same as 
those in Fig Q] except for h = 1.0. One can see that as a 
increases R reaches a maximum R c at a ~ 0.031 and then 
decreases. Obviously, there exists a maximal nucleation 
rate that occurs at a moderate degree of network modu- 
larity. In Fig[3fb), we plot the results of the nucleation 
rates, i?i and i?2, for two-steps process as functions of 
a. As a increases, R\ seems to exponentially decreases 
with a, while i?2 increases monotonically until a = 0.051 
is reached. For a > 0.051, nucleation becomes one-step 
process so that i?2 can not be well defined and the overall 
nucleation rate is only determined by R\. From Figj3jb) 
one can find that R2 is much lower than R\ when the 
value of a is relatively small, so that R is dominantly 
determined by the second step nucleation. While for 
a > 0.031, R is almost determined only by the first step 
nucleation. Thus, there exists a region 0.001 < a < 0.031 
where R is determined by both i?i and R2 ■ Note that we 
have also made extensive simulations for other param- 
eters such as h — 0.7,1.2 and T — 1.5,1.8, and found 
that the qualitative features of the above results do not 
change (results now shown). 

To further understand the above results, we calculate 
free energy of the system by US method, in which we have 
used a bias potential 0.1fceT(A — A) 2 , with A being the 
center of each window. The free energy AF as a function 
of A for three different values of a are depicted in FigH^a) . 
For a — 0.001 and a = 0.031, there are two free-energy 
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FIG. 3: (Color online) (a) The logarithm of nucleation rare 
InR as a function of the degree of modularity a. (b) InRi 
(squares) and InRi (circles) as a function of cr, and dotted 
line indicates the overall rate InR. The parameters are same 
as in FigUJexcept for h = 1.0. 




50 100 150 200 250 300 350 

X 

95 i i i i i i i i i i i i i i i i i i i i i 

a 



90 
85 
80 



■(b) 



j i L. 



_i i i_ 



260 

250*^ 

240 



0.00 0.02 0.04 0.06 0.08 0.10 




0.00 0.02 0.04 0.06 0.08 0.10 
a 



maximums, occurring at the locations of critical nucleus 
of A = X\ and A = AJ, respectively. This picture is 
consistent with the two-step nucleation process described 
before. For a larger a = 0.101, just the first free-energy 
barrier is present, implying that the nucleation becomes 
one-step process. With the increment of cr, A| moves to a 
larger value while the value of A 2 gets smaller, as shown in 
FigHJb). FigBJc) shows that the first free-energy barrier 
Ai*\*, defined as the difference between the free energy 
at A| and the first minimum in free energy (A = 23), 
nearly increases linearly with cr, while the second free- 
energy barrier AF 2 * (definition is similar to that of AF£ , 
and the second minimum in free energy is an increasing 
function of cr, within the range A € [78,93]) decreases 
monotonically with a until AF 2 * vanishes at cr > 0.051, 
which is in agreement with the result of Figj3jb). 



IV. MEAN FIELD ANALYSIS 

In order to unveil possible mechanism behind the above 
phenomenon, we will present an analytical understand- 
ing by CNT and simple MF approximation. Firstly, let 
us assume, for the first-step nucleation, that A nodes 
are in up-spins and the remaining nodes are in down- 
spins in one of modules (say module I for convenience), 



FIG. 4: (Color online) (a) Free energy AF as a function of 
A for three different a = 0.001,0.031,0.101. For smaller a 
two free-energy barriers are clearly observed, while for larger 
a the second one vanishes, (b) The size of critical nucleus Ai 
and A2, and (c) the free-energy barriers AFf and A_F 2 *, as 
functions of cr. The other parameters are the same as Fig[2] 



and all the nodes in the other module (module II) are 
in down-spins. The energy change due to the spin-flip 
of these A nodes can be expressed as the sum of two 
parts AUi = —2h\ + 2JN\ n , where the first part de- 
notes the energy loss due to the creation of A up-spins, 
which favors the growth of the nucleus, while the sec- 
ond part denotes the energy gain due to the forma- 
tion of N{ n new interfacial links between up and down 
spins, which disfavors the growth of the nucleus. Ac- 
cording to MF approximation, N{ n can be written as 
N\ n — pi A ( y — A) + poA-^, where the first part and the 
second part arise from interfacial links inside module I 
and between modules, respectively. For the second-step 
nucleation, we assume that all the nodes in module I are 
in up-spins, and A nodes are in up-spins and the remain- 
ing nodes are in down-spins in module II. This process 
creates new interfacial links inside module II, and at the 
same time removes old interfacial links between module I 
and module II. Thus, the energy change for this process is 
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FIG. 5: The results of mean field analysis, (a) The size of 
critical nucleus A* and A2, and (b) the free-energy barriers 
AFi and AF£ , as functions of a. The other parameters are 
the same as Fig[2] 
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j 2 = -2h\ + 2JN* n where N l 2 n = Pl X (f - A) -p D Af 
is the net number of the interfacial links. The en- 
tropy changes for the two nucleation processes are both 
AS = -h^L [f In (§) + (1 - f ) In (l - § )] . Then, 
the changes of free energy for the two-step processes are 
AF t = AUi - TAS (i = 1,2). In Fig[5]we give the ana- 
lytical results of the critical nucleus and free-energy bar- 
riers as functions of the network modularity. Clearly, the 
analysis qualitatively agrees with the simulation results 



of Fig|4j From Fig|5j one can see that with the increment 
of (7 the size of the first critical nucleus and the height of 
the first free-energy barrier increase almost linearly, while 
the size of the second critical nucleus and the height of 
the second free-energy barrier decrease until a ~ 0.13 is 
reached. This implies that the analysis also predicts the 
extinction of the second nucleation stage, but this pre- 
diction obviously overestimates the transition value of a. 



V. CONCLUSIONS 

In conclusion, we have studied nucleation dynamics of 
Ising model in modular networks consisted of two ran- 
dom networks. Using FFS method, we found that as 
the network modularity gradually becomes worse a tran- 
sition occurs from one-step to two-step nucleation pro- 
cess. Interestingly, the nucleation rate is a nonmono- 
tonic function of the degree of modularity and a max- 
imal rate exists for an intermediate level of modularity. 
Using US method, we obtained free energy profiles at dif- 
ferent network modularity, from which one can see that 
two free-energy barriers exist for very good modularity 
and the second one vanishes when the network modu- 
larity becomes worse. This picture further confirms the 
FFS results. Finally, a mean field analysis is employed 
to understand the nature of nucleation in modular net- 
works and the simulation results. Since stochastic fluctu- 
ation and the coexistent of multi-states are ubiquitous in 
social and biological systems, our study may shed valu- 
able insights into fluctuation-driven transition phenom- 
ena taking place in network-organized systems with mod- 
ular structures. 
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